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We describe the relation between the isothermal atomic compressibility and density fluctuations 
in mixtures of two-component fermions with population or mass imbalance. We derive a generalized 
version of the fluctuation-dissipation theorem which is valid for both balanced and imbalanced 
Fermi-Fermi mixtures. Furthermore, we show that critical exponents for the compressibility can be 
extracted via an analysis of the density fluctuations when phase transitions occur as a function of 
population imbalance or interaction parameter. In addition, we demonstrate that in the presence 
of trapping potentials, the local compressibility, local spin-susceptibility and local density-density 
correlations can be extracted from experimental data via a generalized local fluctuation-dissipation 
theorem which is valid beyond the local density approximation. Lastly, we use the local density 
approximation to calculate theoretically the local compressibility, local spin-susceptibility and local 
density-density fluctuations to compare with experimental results as they become available. 

PACS numbers: 03.75.Ss, 03.75.Hh, 05.30.Fk 



X 



I. INTRODUCTION 

Very recently experimental advances in Bose and Fermi 
systems have allowed for studies of density fluctuations 
and the use of the fluctuation dissipation theorem to 
obtain information about some thermodynamic prop- 
erties of ultra-cold atoms. In the fermion case, the 
measurement of density fluctuations and of the atomic 
compressibility was performed for non-interacting three- 
dimensional systems in harmonic traps [11, , while in the 
boson case, the connection between density fluctuations 
and compressibility was used to study superfluidity in a 
two-dimensional system, and extract critical exponents 
associated with the transition from a normal Bose gas 
to a Berezinskii-Kosterlitz-Thouless superfluid 0. The 
experimental extraction of the isothermal compressibil- 
ity from measurements of density fluctuations was sug- 
gested several years ago both in harmonically conflned 
systems Q and optical lattices [sj , but only recently im- 
provements in the detection schemes of density fluctua- 
tions became sufficiently sensitive to extract this infor- 
mation from experimental data [6l4l0j. 

We see no major technical impediment to use tech- 
niques that are sensitive to density fluctuations in pop- 
ulation imbalanced Fermi-Fermi mixtures with equal 
masses or with unequal masses where, in 

principle, the compressibility and spin susceptibility ma- 
trix elements can be directly extracted from the density 
and density fluctuation proflles. In fact, in a very recent 
experiment using laser speckles, the isothermal com- 
pressibility and the spin susceptibility were measured as 
a function of interaction parameter via the fluctuation- 
dissipation theorem throughout the evolution from BCS 
to BEC superfluidity in balanced Fermi systems. 

In this paper, we extend some of our recent results (l6j . 
and show that through local measurements of densities 
and density fluctuations several local and global thermo- 
dynamic properties can be extracted. We derive a gener- 



alized version of the fluctuation-dissipation theorem valid 
for any mixtures of atoms, and use it to analyze density 
fluctuations and the compressibility of mixtures of two- 
component (or two types of) fermions with and without 
population imbalance at low temperatures. For spatially 
uniform systems, we show that the global compressibil- 
ity can be extracted from measurements of the density 
and density fluctuations for each component, while for 
spatially non-uniform systems, we show that the local 
and global compressibility can be extracted from mea- 
surements of the local density and local density fluctu- 
ations for each component. Finally, we also show that 
across the boundaries between different phases critical 
exponents can be extracted provided that there is suf- 
ficiently spatial resolution for the critical region to be 
experimentally accessed. 

This paper is organized as follows. In section |lTl we 
describe the hamiltonian corresponding to balanced or 
imbalanced Fermi mixtures of equal or unequal masses. 
In section Hm we analyze the compressibility matrix, and 
derive a generalized fluctuation dissipation theorem relat- 
ing the compressibility matrix and density fluctuations 
for balanced or imbalanced two-component Fermi sys- 
tem. Furthermore, we establish a connection between 
the isothermal compressibility or the spin susceptibility 
to the compressibility matrix elements. In section ITVl we 
discuss the effects of traps and derive a local generalized 
ffuctuation-dissipation theorem which is valid beyond the 
local density approximation (LDA). Using this theorem, 
we describe the local isothermal compressibility and the 
local spin-susceptibility, and their relation to the local 
compressibility matrix. In Section|Vl we apply the results 
obtained in the previous sections to describe Fermi- Fermi 
mixtures of equal or unequal masses, as well as equal or 
unequal populations without a trap. In particular, the 
compressibility and spin-susceptibility are calculated as 
a function of population imbalance and interaction pa- 
rameter for the case of mixtures of equal mass fermions. 
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In Section IVIl we use the local fluctuation-dissipation 
theorem and the local density approximation to calcu- 
late the spatial dependence of the local compressibility 
matrix, local isothermal compressibility and local spin- 
susceptibility for population balanced and imbalanced 
mixtures of fermions with equal masses. In section [Vll\ 
wc summarize our main findings and provide some con- 
cluding remarks. 



II. HAMILTONIAN 

To investigate the physics described above, wc start 
with the real space Hamiltonian (?i = 1) 



H ^ Hq + Hfca 



(1) 



for a three dimensional s-wave superfluid. Here, Hq de- 
scribes the kinetic energy of fermions 



2m„ 



Mr), (2) 



where V'i (r) represents the creation of fermions of species 
a with mass , and is the chemical potential which 
determines the average number of fermions for each 
species a. The term that contains the trapping potential 
is 



Htrap = ^ j drVa (r) ha (r) , 



(3) 



where nc,(r) = V'^(r)'0c(r) is the particle density opera- 
tor of fermions of species a, and Va(v) is the correspond- 
ing trapping potential which may vary from one species 
to another. Since we are confining ourselves to s-wave 
interactions the interaction term is 

Hint = / rfrrfr' V^nt(r,r')n^(r)n;(r'), 



involving only spin-up and spin-down densities. In addi- 
tion, we consider the case of an attractive point contact 
interaction T4nt(r,r') = —gS{r — r') with g > 0. 

By allowing each fermion type to have their own chemi- 
cal potential and mass ma , we can investigate thermo- 
dynamic properties of population imbalanced mixtures 
of Fermi systems with equal masses to-|- = mi = m, and 
more generally, unequal masses m^ ^ m^. From the 
grand partition function 



(4) 



Z = Trexp [-H/Tj , 
we can define the thermodynamic potential 
n = -T\nZ, 



(5) 



and obtain global thermodynamic properties such as the 
average number of particles, entropy, specific heat, com- 
pressibility, and so on. In Fermi-Fermi mixtures, the 



compressibility is an important global thermodynamic 
property which can now be directly measured experimen- 
tally in the context of ultra-cold atoms, thus we discuss 
it next. 



III. COMPRESSIBILITY MATRIX WITHOUT 
TRAPPING POTENTIAL 

In this section, we discuss the compressibility matrix, 
the isothermal compressibility, and the spin susceptibil- 
ity for population balanced and imbalanced Fermi-Fermi 
mixtures in two situations. First, we ignore the trapping 
potential Va(r), and analyze the spatially homogeneous 
system to simplify the discussion. After that, we add 
back the trapping potential Va (r) , and return to the spa- 
tially inhomogeneous system, which is more relevant to 
experiments, as described in Sec. IIVI 



A. Compressibility Matrix 

Setting the trapping potential Va (r) to zero in Eq. ^ 
leads to no-trap Hamiltonian 



7? = 7?i -^/i„iV„, 
where the first term is given by 



(6) 



V2 



2ma 



Mr) + Hint, (7) 



representing the sum of the kinetic energy and interaction 
contributions, and Na = / drna{v) the number operator 
for hyperfine state a. With this explicit separation, it 
becomes easier to see that the average number of particles 
Na = (Na), defined by the thermodynamic average 



(Na) = Z-^Tr 



Nae 



(8) 



can be rewritten in terms of the thermodynamic potential 
51 as 



Na 



{Na 



dHa / rp 



(9) 



Starting from Eq. we define the pseudo- 

compressibility matrix {k} through its matrix elements 



rrfdNa\ 



= -T 



\d^ladflp 



(10) 



Using the fact that Na commutes with Hi defined in 
Eq. ([7]) and using the definition of in Eq. ([5]), we can 
calculate the matrix elements of {k} explicitly, leading 
to 



nc.p^{NaNp)-{Na){Np), 



(11) 
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expressed in terms of thermodynamic averages. Here, the 
correlated average of number operators and Np is 



(N^Np) = Z-'Tt \Na.Np e-"/^ 



(12) 



The mechanical stability of the system at any fi- 
nite temperature T is guaranteed when both eigen- 
values of the matrix {k} are positive definite, while 
in the limit of T — > it is the positivity of both 
eigenvalues of \imT^o{{'^} /T) that guarantees the sta- 
bility of the system. At finite T, the eigenvalues of 
{k} are positive if the determinant det[{K}] > and 
the trace Tr[{K}] > 0. Analogously, at T ^ 0, 
the eigenvalues of \imT^o{{K.} /T) are positive if the 
determinant det [limT^o({'«}/r) ] > and the trace 
Tr[limT^o(W/T)] >0. 

Using the definition Na = {No) , we can explicitly iden- 
tify that the pseudo-compressibility matrix is a measure 
of the correlation in density-density (or number-number) 
fluctuations, since it can be rewritten as 



(13) 



where ANa = Na — Na- The corresponding generalized 
compressibility matrix {k} can be obtained from {k} 
through the relation 



«Q/3 = liap/iNaNp), 



(14) 



describing a generalized fluctuation-dissipation theorem 
for multicomponent fermions, which is valid for fermions 
of equal or unequal masses and equal or unequal popu- 
lations. As it is discussed below, when the populations 
for each fermion species are the same Nf = N^, this 
more general result reduces to the standard form of the 
fluctuation-dissipation theorem [l7| . 

In order to be more explicit, we note that the com- 
pressibility matrix element is directly related to the 
particle number fluctuation in the spin- up/spin- up chan- 
nel ((AiV-t-)^), that K^i is directly related to the particle 
number fluctuation in the spin-down/spin-down channel 
{{ANi)"^), and that k^^ is directly related to the cor- 
relation between the particle number fluctuation in the 
spin- up/ spin- down channel (ATV^AiVj,). In addition, the 
matrix {k} is symmetric, such that /t-i-j. = k^-^. 

In very recent experiments [HQ, particle number (den- 
sity) fluctuations were directly observed for an ideal (non- 
interacting) Fermi gas of equal masses and equal popula- 
tions = N^, in which case the chemical potentials fia 
for each hyperfine state a are exactly the same {fif = /ij.). 
For arbitrary interactions, identical masses and popula- 
tions, it is sufficient to analyze fluctuations only in the 
total average number of particles N = N^ + Ni to extract 
the isothermal compressibility 



1 fdV 



V \dP 



V_ fdN_ 



1^] =T72 hr- • (15) 



Using the standard form of the fluctuation-dissipation 
theorem for a balanced system, the isothermal compress- 
ibility can be rewritten as 




(16) 



For mass or population imbalanced systems, how- 
ever, the relation between the isothermal compressibility 
kt = —{dV/dP)T/V and the fluctuations in the number 
of particles is more subtle. The determination of kt h> 
volves an analysis of the pseudo-compressibility matrix 
{k}, which contains all the information about fluctua- 
tions in particle numbers (densities) for each hyperfine 
state. Using similar experimental techniques as those 
described in Refs. [J 0], it should be possible to mea- 
sure the matrix elements of {/t} directly not only for 
the ideal non-interacting Fermi system, but also when 
strong correlations are present, as the fluctuation dissi- 
pation theorem described above is valid for all interaction 
strengths. Thus, it is important to identify the relation 
between the isothermal thermodynamic compressibility 
Kt and the elements of the compressibility matrix {k}. 
This connection is discussed next. 



B. Isothermal compressibility 

The relation between kt and Kq/j can be established 
via the thermodynamic potential £7 = —PV, where P is 
the pressure and V is the volume of the system. Defining 
G = fl + PV = 0, and recalling that is a function of 
temperature T, volume V and chemical potentials fj,-f, fi^ 
results in 

dG = -SdT + VdP - Npdnfs = 0. (17) 

At constant temperature dT — 0, we can establish the 
relation 



(18) 



This means that the inverse isothermal compressibility 
(also called the Bulk modulus) 



1 



dP_ 
dV 



(19) 



can be written in terms of isothermal partial derivatives 
oi jip with respect to volume 



(20) 



Ni. 



But, in turn, the partial derivatives {dfj.p/dV)T can be 
expressed in the terms of the average numbers N-^ and 
Ni, since /i^g = ^^{N^/V, Ni/V). For the isothermal 
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variation of the spin-/? chemical potential with respect to 
the volume V, we can write 



(21) 



The last step in establishing the relationship between 
Kt and Kq^, just requires the use of the definition 
{dNaldni3)T = Kq^/T from Eq. (fTO|) or its inverse 
{d^ip / dNa)T = T/kap. Putting all this together leads 
to the interesting relation 




(22) 

which can still be further written in a compact form 



afj 



where we used the definition of the compressibility ma- 
trix elements Kai3 = Kap/iNaNjs) described in Eq. ([H]). 
From the relation between the pseudo-compressibility 
matrix elements kap and the density-density (number- 
number) fluctuations shown in Eq. (|13p . we can rewrite 
the inverse of the isothermal compressibility as 



- = -E 



(24) 



The result given in Eq. reminds us of the for- 
mula for the computation of the equivalent resistance of 
four different resistors in parallel. Each matrix element, 
1/kq;3, acts as independent parallel channel of number 
fluctuations, which combined give essentially the inverse 
isothermal compressibility. Measurements of the com- 
pressibility matrix elements are important as they 
can be connected further to the pseudo-spin susceptibil- 
ity of the system, as discussed next. 



C. Connection to pseudo-spin susceptibility 

At this point it is important to mention that even 
though we have chosen to work with particle numbers 
iV-|-, and chemical potentials fi-^, fi^^ we could have 
chosen also to work with the total number of particles 



N+=N^ + N^, 
the difference in particle numbers 



(25) 



(26) 



and their corresponding chemical potentials /i± = (/i-|- ± 
/ij^)/2, respectively. This new choice provides the same 
information about the system since the term involving 



the chemical potentials in the Hamiltonian H of Eq. 
can be rewritten as 



(27) 



such that the corresponding partition function Z or ther- 
modynamic potential also contain the same terms. In 
this case, we can define a similar pseudo-compressibility 
matrix 



T 



-T 



_dH}_ 



(28) 



where each of the indices i,j can take ± values. An ex- 
plicit calculation using the grand partition function Z 
allows to express kij in terms of the thermodynamic av- 
erages 



(29) 



The connection between the pseudo-compressibility ma- 
trices in the two different representations is simple. The 
first diagonal term is = k^-f + k^ + 2k-fi, the sec- 
ond diagonal term is k = k-f^ + k^ — 2k-^^, while the 

off-diagonal terms k-| = k |_ = k-^-^ — k^^ are identical 

since the matrix kij is symmetric. 

The corresponding expression of the isothermal com- 
pressibility Kt has exactly the same form as before, just 
with indices (t,i) mapped into indices (+, — ), and thus 
leading to 



1 

Kt 



T 
V 



k++ 



k+- 



k-+ 



7V2 



(30) 



This way of writing kt is more transparent because 
the population balanced case iV-t- = Ni corresponds to 
N- ~ 0, and thus leads to the standard expression for 
the compressibility 



1 



T f Nl 



V \ k 



(31) 



which upon inversion and further identification of A'^-i- — >■ 
N leads to the balanced form of the fluctuation dissipa- 
tion theorem as described in Eq. (fT6|) . 

For mass or population imbalanced cases, the direct 
use of Eq. ([50)) allows for the extraction of the atomic 
compressibility directly from measurements of densities 
and density-density fluctuations. This result general- 
izes initial suggestions that the isothermal compressibil- 
ity could be measured experimentally for balanced sys- 
tems in optical lattices or harmonic traps 0, [E] via the 
fluctuation dissipation theorem. 

Notice that the dimensionless pseudo-spin susceptibil- 
ity Xzz of the system can also be extracted from the 
pseudo-compressibility matrix kij^ since (/i-j- — 

/i^)/2 plays the role of an effective magnetic fleld hz along 
the quantization axis z, and N- = N^ — N\^ plays the role 
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of the magnetization m^. For imbalanced Fermi systems, 
this implies that 



T 



fdN_ 



T 



( dm. 



Tx.z. (32) 



and gcnerahzcs the results obtained for ultra-cold 
fermions interacting via p-wave interactions (triplet chan- 
nel) [13. 

Now that the case of zero trapping potential has been 
analyzed, we add back Va{Y) to the Hamiltonian of the 
system and discuss its effects next. 



IV. COMPRESSIBILITY MATRIX IN A TRAP 

To consider the most general case, where the external 
trapping potential V^(r) is non-zero, we include a source 
term in the Hamiltonian defined in Eq. (IH), leading to 



(33) 



In the presence of this source term, the grand partition 
function Z becomes a generating functional of the Jq, 
as well as temperature T, volume V , and the chemical 
potential jia, 

Z[T,V, Ma, Ja(r) ] = Tie-"'/^. (34) 

The thermodynamic potential 

= -TlnZ[T,F,Aia,Ja(r)] 

is also a functional of Jq (r) , such that local and even non- 
local thermodynamic quantities can be obtained. For in- 
stance, the local density ^^(r) = {na(r)) is simply writ- 
ten as the functional derivative 



na(r) = -T 



sn 



5 Mr) 



(35) 



And the second derivative gives the correlation of the 
particle density fluctuation 



(An„(r)An^(r')) =-r 



,52 r! 



5Ja(r) 5Jp{v') 



(36) 



with Ana{r) = ^^(r) — na{r). Correspondingly, we de- 
fine the non-local pseudo-compressibility matrix elements 



KQ/3(r, r') = {Aha{r) An^(r')). 



(37) 



Then, the local pseudo-compressibility matrix is simply 
defined as kajsir) = kai3{r,r). 

The inverse of the local isothermal compressibility 
KT(r) can be expressed in terms of the local pressure 
density p{r) as 



1 



KT(r) 



= -V 



dp{r) 
dV 



(38) 



where the total pressure is P = / drp(r). Following the 
same steps used to derive Eq. ([^5]) , we arrive to the local 
relation 



1 



T 



'«T(r) V ^ Ka/3(r) 
where the local compressibility matrix becomes 

'«a,9(r) 



KQ/3(r) 



(39) 



(40) 



thus revealing a local generalization of the fluctuation- 
dissipation theorem. 

Analogously to the zero trap case, instead of working 
with the local particle number densities nQ,(r), we could 
have chosen to work with the total number of particles 

n+(r) = 7it(r)+?i^(r), (41) 

and the difference in the number of particles 

n_(r) = nt(r) -n;(r). (42) 

In this new basis, we define a similar pseudo- 
compressibility matrix 



K,y(r,r') = (An,(r) Anj(r')). 



(43) 



The corresponding expression for the local isothermal 
compressibility KT(r) has exactly the same form as be- 
fore, just with indices (t,!) mapped into indices (+, — ), 
thus leading to 



1 _ ^ \p 1 
where the local compressibility matrix elements are 



Ky(r) 



71, (r) nj(r)' 



(44) 



(45) 



The use of equation Eq. (|44p allows for the extraction of 
the local isothermal compressibility KT(r) and the local 
dimensionless pseudo-spin susceptibility 



^__(r) _ ((An_(r))2) 



T 



T 



(46) 



directly from the measurements of the local densities and 
local density fluctuations for imbalanced systems. 

Now that the case of non-zero trapping potential has 
been analyzed, we turn our attention to the application 
of these results to the cases of balanced and imbalanced 
Fermi-Fermi mixtures. 



V. FERMI-FERMI MIXTURES WITHOUT A 
TRAP 

So long as the particle number operators commute 
with the Hamiltonian of the system, the results discussed 
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in section IIIII apply to any binary mixtures (Bose-Bose, 
Fermi-Bose, Fermi-Fermi) in the case where there is no 
trap; while the results of section ITVl apply to any binary 
mixtures in a trap. In this section, we discuss generally 
the cases of Fermi-Fermi mixtures of equal or unequal 
masses when populations are balanced or imbalanced, 
but we focus particularly in the case of equal mass mix- 
tures with balanced or imbalanced populations as a con- 
crete example. 



A. Thermodynamic Potential without a Trap 

As an specific example of the general relations just 
derived, wc discuss the case of imbalanced Fermi-Fermi 
mixtures with equal masses, which has attracted substan- 
tial interest p9l - l24| recently. 

The action corresponding to the Hamiltonian given in 
Eq. ([1]) is simply written in the coherent state represen- 
tation as 

S = j^^ dr j dr [i^^{T, T)dMr, r) + H(r, r)] , (47) 

where '0Q,(r,T), and ?/'c((r,T) are Grassman variables for 
the fermion type a. The corresponding partition function 
is then written as the functional integral 



a. 

] exp ( — 5 [V'a, V'a] ) 



(48) 



from which the thermodynamic potential fi = — T In Z 
can be directly calculated. 

Since we are considering first the case without a trap- 
ping potential, Va{v) = 0, we can use the transla- 
tional invariance of the system to write the Hamiltonian 
H{t) = / rfr-H(r,r) as 

H{r) = ^Ck,cVik,a(T) Vk,a(r) -g5^Bq(r)Sq(r) 

k,a q 

(49) 

where ^k,a = ek,a — ^J■a, with ek,Q = k'^/{2ma) being the 
single particle dispersion. Here, ipk^aiT) is the Fourier 
transform of '(/'^(r, r), and 

Bciir) = XI V^q/2-k4(^) V'q/2+k,t(T-) (50) 
k 

describes paired fermions and Fk corresponds to the sym- 
metry of the pairing interaction. For instance, in the case 
of s-wave pairing we have Fk = 1. 

The effective action S has been successfully calculated 
for the case of a uniform superfluid f25l| in the Gaussian 
approximation as 



the inverse fluctuation propagator. The explicit expres- 
sion of F~^(q) is given in Appendix [K\ 
The saddle point action 5*0 is given by 



So 



^X(a,+ -a,-)+ J2 inK(-£:k,.)], 



where fk.i. = ^k,- + (-l)'' ( Ck,+ + l^kP are the quasi- 
hole and quasi-particle energy spectrum for v = I and 
1/ = 2, respectively. In addition, Ak = AgFk is the order 
parameter for superfluidity for pairing with zero center of 
mass momentum, np{E) = 1/(1 + exp ii'/T) is the Fermi 
distribution, and ^k,± = (^k,t =t ^k.i) /2 = fc^/2m± — /i±, 
where m± = 2m^mi/ (mj, ± m-|-) and /i± = (/i-|- ± //j,) /2. 
Notice that m+ is twice the reduced mass of the t and 
\, fermions, and that the equal mass case (m-|- = m^) 
corresponds to |m_| — >■ oo. 

The fluctuation term in the action leads to a correction 
to the thermodynamic potential, which can be written as 



Via — + r^fiiict, 



(51) 



where flo = TSq and ^f^^^t = TJ2q Indet [TF-^{q)] . 

The saddle point condition SSq/SAq = leads to the 
order parameter equation 



E 



I I T,' — 



(52) 



where - {^1^ + \A^\^)^^^ and = n;^(£k,i) ± 
np[£\^^2)- As usual, we eliminate g in favor of the scat- 
tering length via the relation 



1 
9 



47ra„ 



E 



iFkl^ 



ek,t + £k,i 



The order parameter equation needs to be solved 
self-consistently with the number equations = 
— {dfl/d^a)T, which has two contributions 



(53) 



Here, No^a = —{d^a/d^a)T is the saddle point number 
equation given by 



E 

k 



— -^k 



(54) 



where Sq. = 7q(1 — X^) with 7-f- = 1 and 74^ = —1. And 
A^fluct,a = — (i917fiuct/(?/ia)|T IS the fluctuation contribu- 
tion to A'q, given by 



d 



.a = -rV — Indet 



TF-\q) (55) 



where q = (q, t-f), with = 2TTiT being the bosonic Mat- 
subara frequency at temperature T. Here, A{q) is the or- 
der parameter fluctuation field and the matrix F~^{q) is 



In order to extract the isothermal compressibility kx 



of the system or the bulk modulus B 



shown in 



Eq. (j22p and the pseudo-spin susceptibility x^z defined 
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in Eq. (|32)) , we need to obtain the pseudo-compressibility 
matrix {k} described in Eqs. ([TU|) and ([T^ . Using 
the self-consistent solutions from the order parameter 
Eq. ([5^ and the number equation, Eq. ([55)1 . leads to a 
thermodynamic potential f2 which is a function of the 
order parameter Aq, the chemical potentials fia, and 
temperature T. However, from a thermodynamic point 
of view fi is ultimately dependent on fia and T only, 
such that [Ao(Ma,T),^„(r),r] ^ [/Zc(T), T] . This 
implicit dependence of the order parameter Aq on the 
chemical potentials fi^, fi^, and temperature T requires 
that derivatives of the thermodynamic potential n in- 
clude both their explicit and implicit dependencies on 
chemical potentials, and temperature. This is extremely 
important for the calculation of the isothermal compress- 
ibility, spin susceptibility and pseudo-compressibility ma- 
trix to be discussed next. 



B. Pseudo-compressibility matrix without a trap 

Because of the implicit dependence of Aq on ^q,, the 
direct calculation of the pseudo-compressibility requires 



N 


/ NU / 


U 


1 V 


\ NTI \ 


TT 
U 


I -1 


1 2 


3 ^ 




l/{kF+as) 





0.5 



^ 



-0.5 



-1 



FIG. 1: The zero temperature phase diagram of popula- 
tion imbalance P = versus scattering parameter 
l/{kF+as) for equal mass fermions is shown. For P 7^ 0, 
the normal (N) state becomes superfluid (denoted by NU 
and U) as the scattering parameter l/(fcF+as) increases. U 
represents the uniform superfluid phase, and NU represents 
a non-uniform superfluid region, where the U phase is unsta- 
ble. For P — 0, the U phase is stable for all values of the 
scattering parameter. 



Kg/a 

T 



T,c 



9|Ao| 



T,c 



9|Aof 



(56) 



T,i 



where the label "e" ("i") means explicit (implicit) deriva- 
tive. The explicit expressions of the matrix elements at 
zero and finite temperature are given in Appendix |B] For 
any non-zero temperature T, the mechanical stability of 
the uniform superfluid and normal phases is guaranteed 
if all eigenvalues of kap are positive, or, equivalently, if 
all eingenvalues of kij are positive. Since we are dealing 
with 2x2 matrices, this is achieved when Tr [{k}] > 
and det [{k}] > 0. When the lowest eigenvalue of {k} 
reaches zero then the system becomes mechanically un- 
stable. The existence of mechanical stability coincides 
with the condition that the thermodynamic potential is 
a minimum [Bj with respect to variations of the order pa- 
rameter |Aop, this means that it satisfies simultaneously 
the extrcmum condition (9r2/9|Aop) = 0, and the posi- 
tive curvature condition (9^r2/9|Ao|^) > 0. For zero tem- 
perature (T = 0), the mechanical stability is guaranteed 
when both the eigenvalues of the matrix Iuht^o [{'^l/^^l 
are positive. 

Using the stability conditions described above in com- 
bination with the solutions of the order parameter 
Eq. (|52p and number Eq. (|53p , we construct the zero tem- 
perature phase diagram for a mixture of fermions of equal 
masses, different hyperfine states and no trapping poten- 
tial. The zero temperature phase diagram of population 
imbalance P = {N^ - 7V;)/(iV-^ + N^) = N^/N+ versus 
scattering parameter 1/(A;f+Os) is shown in Fig.[TJ Here, 
kp_^_ = {kp^ + kp^)/2 is an effective Fermi momentum, 
where kpa is the Fermi momentum of each species, and 
eF+ = kp_^/ (2m+) is an effective Fermi energy, which 
fixes the energy scale. The effective Fermi momentum 



kp^ is convenient because it fixes the total number of 
particles = + N^, where the density for particle 
a is {N^/V) - klJ{6TT^). 
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FIG. 2: The pseudo-compressibility matrix elements kafs/T 
at T = are shown in a) as a function of population imbalance 
P for interaction parameter l/(fcF+as) = 1.92, and in b), as a 
function of l/(fc_F+as) for P = 0.5. The matrix elements in the 
region where the uniform superfluid phase is unstable is shown 
in gray. The inset in b) shows an expected discontinuity at the 
boundary between the normal and the non-uniform superfluid 
state. Notice that, in the normal state, k^i = k^^ = 0. 

In Fig. [21 we show the pseudo-compressibility ma- 
trix elements kap in a) as a function of population im- 
balance P for fixed l/{kp^as) = 1.92 and in b) as 
a function of l/{kp^as) for fixed population imbalance 
P = 0.5. As seen in Fig. when population imbalance 
is changed in the BEC regime, e.g. l/{kp^as) = 1.92, all 
limT^o i^ap /T diverge at the phase boundary between 
the uniform superfiuid phase and the non-uniform phase 
as kaf) - (P-Pc)-" with P^ = 0.5 for l/ikp+as) = 1.92. 
The critical exponent associated with this divergence is 
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ly = 1 and can be extracted from the density fluctua- 
tions revealed in the matrix elements of Kq.^. In Fig. ^jp, 
the matrix elements limT^o i^ap/T arc shown as a func- 
tion of the interaction parameter l/kp^as for fixed pop- 
ulation imbalance P = 0.5. As shown in the inset of 
Fig [2b, discontinuities are expected when the interaction 
parameter l/{kp^as) is increased to cross the boundary 
between the normal phase and the region where a non- 
uniform supcrfluid phase is present. The negative values 
for the matrix elements lim7^_>o i^ap /T just indicate a re- 
gion of non-uniform superfluidity, that is, a region where 
uniform superfluidity is not mechanically stable. When 
the boundary between the non-uniform phase and the 
uniform superfluid phase is crossed then limT-s-o K,ap/T 
diverges as (A — Ac)~^, where A ~ l/^kp^cis), and Ac is 
the critical interaction parameter. For population im- 
balance P = 0.5, the critical interaction parameter is 
Ac = 1.92. This illustrates the point that using the 
generalized fluctuation-dissipation theorem described in 
Eqs. and allows for the extraction of critical 
exponents of density-density correlations accross phase 
boundaries. Notice that in this approximation and within 
the normal state region, all limr-^o IT arc constant 
as a function of interaction parameter Xj (kpj^as)- Addi- 
tional corrections to the present approximation are nec- 
essary to capture fully interaction effects in the normal 
state. 



C. Isothermal Compressibility and Spin 
Susceptibility 

The isothermal compressibility per unit volume kt = 
ktIV and the pseudo-spin susceptibility Xzz are shown 
in Fig. [3^ as a function of population imbalance P for 
= 1-92 as in Fig. [51 Notice that as P in- 
creases, both kt and Xzz diverge as (P — Pc)"" at the 
critical population imbalance Pc = 0.5, and become neg- 
ative for \P\ > Pc signaling a quantum phase transition 
from uniform superfluidity with coexistence of excess un- 
bound fermions and paired fermions in the same spatial 
region into a phase separated regime where excess un- 
bound fermions and paired fermions tend to avoid be- 
ing in the same region of space. The critical exponent 
in the region where the isothermal compressibility and 
the pseudo-spin susceptibility are positive has the sad- 
dle point (mean field) value of = 1. A renormalization 
group calculation would be necessary to obtain correc- 
tions to this critical exponent, but we postpone it to a 
later date. 

In Fig. [3)3, wc show kt and Xzz as a function of 
l/{kp^as) for fixed population imbalance P = 0.5. In 
the uniform superfluid phase, where l/{kp^as) > Ac, the 
same divergent behavior occurs as kt ^ {X — Ac)^^ and 
Xzz ~ (A — Ac)""'^, where A = l/{kp^as) and Ac is the 
critical interaction parameter (Ac ~ 1.92 for P = 0.5). 
The critical exponent v = 1 reflects the saddle point ap- 
proximation used to calculate it. Renormalization group 
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FIG. 3: The isothermal compressibility per unit volume V, 
Rt ~ f^r/V as well as the pseudo-spin susceptibility Xzz = 

k /T in the limit of T = are shown in a) as a function of P 

for 1/ (kp+as) — 1.92, and in b) as a function of l/(fc_F+as) for 
P — 0.5. The unit on the vertical axis is chosen as 3A'^+ /4eF+ • 
The gray curves in a) and b) illustrate the instability region of 
the uniform superfluid phase. In c) and d), the ratio Xzz/^t 
is shown as a function of P and l/{kp+as) in a stable uniform 
superfluid state, respectively. 



calculations are necessary to obtain corrections to ly = 1. 

In Fig. [3J:, we show the ratio of the pseudo-spin sus- 
ceptibility Xzz to the isothermal compressibility per unit 
volume Kp 



Xzz 

hp 



= Ni 



2P 



P' 



(57) 



as a function of |P| with a fixed value of 1/ {kp-^-ag) ~ 
1.92, and Fig.[3Ji shows the ratio as a function l/{kpj^as) 
for P = 0.5 in the uniform superfluid phase. In Figs. [5J; 
and[3Jl notice that the ratio Xzz/^p approaches a finite 
value, as indicated by the dotted line, refiecting the same 
power law divergence of Rp and Xzz at the phase bound- 
ary where the uniform (U) supcrfluid phase becomes un- 
stable. 

For completeness, in Fig. |4l we show plots of the zero- 
temperature isothermal compressibility for the popula- 
tion balanced case corresponding to P = 0. In partic- 
ular, we show separately the contributions coming from 
both the explicit and implicit dependencies on the chem- 
ical potential /i in the expression for the isothermal com- 
pressibility 



(58) 



where N = (N) is the total average number of particles, 
and k — T{dN/dfi)T is the 1x1 compressibility matrix 





( k \ 


V 




T 




~ 7V2 
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for the balanced case, that is, the scalar 
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FIG. 4: (color online) {dN/dfi)T as a function of l/{kF+as) 
for T = is plotted by the blue solid line for the population 
balanced case P = 0. The unit on the vertical axis is chosen 
as N/ep- The green (dot-dashed) line represents the contri- 
bution from the explicit derivative [dN / d^)T,c, while the red 
(dashed) line corresponds to the contribution from the im- 
plicit dependence of |Ao| on /i, that is, from {dN/dfi)T i = 
{^N/^\Ao\^)T,.■{^\Aof/^^i)T^. 

As seen in Fig. |4l the second term in Eq. ((59)) . which 
contains the implicit dependence on the chemical po- 
tential /z, is important both in the intermediate region 
— 1 < l/^kpcis) < 1 around unitarity and in the BEC 
regime. Analytical limits in the BCS and BEC regimes 
are also shown and confirm the importance of the second 
term. The derivations are given in Appendix [C] The 
zero-temperature behavior of k/T = {dN/dfi)T in the 
BCS limit is given by 

/dN\ 3N ( ^\ 

having a positive correction with an upward curvature 
as the interaction parameter increases: 7 = [7-\/2 — 4 4- 
3sinh~^(l)] 8/e^ ~ 1.2519. The corresponding result in 
the BEC limit is given by 

IT-^T ' (62) 

a [I Zep kpQs 

which agrees with the result of the compressibility of a 
weakly interacting gas of bosons with mass tub = 2m, 
and scattering parameter qb = 2as [26| . Thus, it is only 
when the implicit dependence of the order parameter Ag 
on the chemical potential fi that the qualitatively correct 



result in the BEC limit is recovered. This can be gen- 
erally understood by the argument that increasingly at- 
tractive interactions tend to make the Fermi system less 
degenerate and thus more compressible. This tendency 
is better illustrated in the BEC regime where molecular 
bosons emerge as new degrees of freedom and the Pauli 
pressure is weakened with increasing interaction parame- 
ter thus leading to a more compressible Fermi superfluid. 

Before we conclude this section it is worth noting that 
the zero-temperature spin susceptibility Xzz is exactly 
zero in the case of balanced populations {P = 0), since 
for a uniform s-wave superfluid at T = all the fermions 
are paired. 

Having discussed the isothermal compressibility and 
spin-susceptibility for s-wave Fermi superfluids without 
a trap, now we turn our attention to the case where the 
trapping potential is non-zero. 

VI. FERMI-FERMI MIXTURES IN A TRAP 

In the case of a non-zero trapping potential VQ.(r), we 
make use of the local density approximation (LDA), and 
obtain the local thermodynamic potential 

n{r)^n[fi^{r)] (63) 

from the thermodynamic potential in the absence of a 
trap, defined in Eq. ([5T|) . via the substitution 

fla ^ f^a{r) = ^^a-Va{r). (64) 

We note, however, that the generalized fluctuation dis- 
sipation theorem described in Sec. I VII can be used 
to extract the local compressibility and the local 
spin-susceptibility directly from experimental results of 
density-density fluctuations, without invoking the local 
density approximation (LDA). 

A. Density profiles and order parameter 

The use of LDA implies that both the order parameter 
Aq and the number of particles Na defined in Eqs. ([52]) 
and ()53p , respectively, become functions of position r via 
the position dependent chemical potentials /^^(r). As a 
result, we have Ao(r) = Aq [/i-|-(r), /i4_(r)] and ^^(r) = 
na[iJ.^{v), fii{r)] , where na{r) — Na{r)/V is the local 
density for fermions of type a. To simplify our discussion, 
we consider isotropic harmonic trapping potentials 

Vcir) = (65) 

where 7q, = mcW^ with uja being the trapping frequency 
of fermion of type a, and r = |r| is the magnitude of the 
position vector r. 

Following the general procedure of Sec. IIVI it is simple 
to show that within LDA the local density is given by 
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FIG. 5: (color online) a) The particle density profiles 
na{x) and b) the Jacobian- weighted '^'Kx'^na{x) normalized 
by A'^+ = (2m+eF+)"^''^/(37r^), are plotted as a function of di- 
mensionless position x — |r|/rTF for fixed interaction param- 
eter l/(fcF+as) = 3.0 at the population imbalances P = 0, 
P = 0.2, and P = 0.6, respectively. The dimensionless or- 
der parameter Ao(2;) normalized by eF+ is represented by the 
black-dot-dashed line, while the densities nf{x) and ni{x) are 
described by the blue-dashed and red-solid lines, respectively. 
Notice that for P = 0, the densities ni^{x) and 714.(3;) coincide 
exactly. 
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X = |r|/rxF 

FIG. 6: (color online) The profile of pseudo-compressibility 
matrix elements Kap{x)/T axe shown in a) as a function 
of dimensionless position x for fixed interaction parameter 
l/ikF+cbs) = 3.0 at population imbalances P = 0, P = 0.2, 
and P = 0.6. And the corresponding Jacobian-weighted pro- 
files A-KX^ K,ap{x) are shown in b). The unit on the vertical 
axis is chosen as 3A'^+/4eF+. The solid-green line denotes 
the tt element of both pseudo-compressibility and Jacobian- 
weighted matrix. The dotted-red line describes ^4- element, 
and the dashed-blue line describes 4-4- element. 



B. Compressibility and spin susceptibility 



where the thermodynamic potential VL is defined in 
Eq. ([?T|) with the replacement Ha — > ^^a{T^), where /iQ(r) 
is the local chemical potential defined in Eq. (|64)) . 



For the equal mass case, we show in Fig. [5]i the spa- 
tial dependence of the particle density profiles 'na{x) and 
the order parameter Aq (x) as a function of dimensionless 
position X = \r\/rTF, where ttf is the Thomas- Fermi 
radius defined through the condition tp+ = 7_|.r^p/2, 
where 7+ =7^ + 74,- These spatial profiles show that 
superfluidity coexists with excess unpaired fermions, e.g. 
spin-up fermions n^(x)^ but the majority of excess un- 
paired fermions are pushed away from the center of 
the trap. This effect is better seen in Fig. [Sh, where 
the Jacobian-weighted particle numbers A-kx^ na[x) are 
shown for interaction parameter 1/ {kp+as) = 3.0 at pop- 
ulation imbalances of P = 0, P = 0.2, and P = 0.6. 

Given our choice of majority spin-up fermions, the den- 
sity of paired fermions npair(r) is controlled by the mi- 
nority spin-down fermions, becoming npair(r) ~ 2n^(r), 
and the density of excess fermions nc(r) is simply the 
difference between the densities of the spin-up and spin- 
down fermions such that no(r) = — n^{r). Notice 



that 



(r) vanishes beyond the critical radius Vc be- 



yond which the order parameter for superfluidity Ao(r) is 
zero. Furthermore, as population imbalance P increases 
the density of pairs npair(r) also decreases concomitantly 
with a shrinkage of the superfluid region determined by 
Tc = |rc|. In contrast, as P increases, the overall density 
of excess fermions no(r) increases. 



Following the steps outlined in Sec. llVi it is straight- 
forward to show that the matrix elements of the pseudo- 
compressibility matrix can be obtained from the relation 



k^pir) = -T 



(67) 



where f2[/i-(-(r), /i^(r)] is the thermodynamic potential 
within LDA. In Fig. the spatial dependence of the 
pseudo-compressibility matrix elements kap{x)/T are 
shown for the equal mass case as a function of dimension- 
less position x for the same parameters used in Fig. [S] 
Notice that for P = all the matrix elements of kap {x) /T 
in the superfluid phase, x < Xc, coincide and vanish at 
X > Xc- For P ^ 0, however, the matrix elements are split 
at a; < Xc, reflecting the coexistence of the superfluidity 
and excess unpaired fermions, and vanish at x > Xc ex- 
cept for k^^(x)/T due to the unpaired spin-up fermions. 
In Fig.[5jD, the Jacobian-weighted pseudo-compressibility 
matrix elements Airx'^kaisix) /T are plotted, showing that 
these matrix elements decrease in the superfluid region 
(x < Xc) as P increases, since excess fermions are trans- 
ferred to the non-superfluid regions (x > Xc). Notice, 
however, that k^^{x)/T in the normal phase beyond Xc 
increases with P for the same reason. Within LDA, each 
matrix element exhibit a significant discontinuous drop 
when crossing the critical radius Xc- 

In Figs. [7] and IH the spatial dependence of the 
isothermal compressibility kt(x) and that of the spin- 
susceptibility Xzz(x) arc shown, respectively, for the same 
parameters as Fig. [SJ As seen in Fig. [71 kt(x) exhibits 



11 



a discontinuity at the position x = Xc where the order 
parameter Ao(a;) and the minority density ?t-j_(x) vanish. 
Since, within LDA, the isothermal compressibihty is 

1 ^ r ^ / na{x)np{x) \ 

the discontinuity becomes evident as '^■4,(3;) — ^ 0, when 
X Xc- For finite population imbalance P, kt{x) in- 
creases monotonically with x in the superfluid region 
x < Xcj and exhibits a discontinuity at Xc- For the pop- 
ulation balanced case P = 0, the isothermal compress- 
ibility kt{x) also increases monotonically with x and ex- 
hibits a discontinuity at the boundary between the super- 
fluid region and vacuum (where there are no fermions). 
These discontinuities at also clearly seen in 

Fig. Wp^ where the Jacobian- weighted isothermal com- 
pressibility profiles Attx^ kt{x) are shown. 

P = P = 0.2 P = 0.6 
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FIG. 7: a) The zero-temperature isothermal compressibility 
profile per unit volome kt{x) and b) the Jacobian- weighted 
profile 4:Tvx^Rt{x) in unit of 3A^+/4eF+ are shown as a func- 
tion of dimensionless position x for fixed interaction param- 
eter 1/ (kp+as) = 3.0 at population imbalances P = 0, 
P = 0.2, and P = 0.6. 

As seen in Fig. [8l for non-zero P, the pseudo-spin sus- 
ceptibility Xzz{x) also exhibits a discontinuity at the po- 
sition Xc, where the order parameter A(a;) and the mi- 
nority density n^^x) vanish. This discontinuity signals 
the phase boundary between the spatial region where su- 
perfluidity coexists with excess fermions, and the spa- 
tial region where only excess fermions exist. In Fig. [5Jd 
the Jacobian- weighted pseudo-spin susceptibility profiles 
Attx'^Xzz{x) are shown. Notice that the overall response 
Xzz = (''tf/^) / dxi'Kx'^Xzzix) increases monotonically 
with P, as more unbound excess fermions are available. 
In the balanced case, where P = 0, we have trivially 
Xzzix) — 0, since for x < Xc all spin-f and spin- J, fermions 
are paired into a singlet state, and for x > Xc there are 
no fermions. 

Now that we have completed our discussion of the spa- 
tial dependence of the particle density, order parame- 
ter, compressibility and spin susceptibility for balanced 
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FIG. 8: a) The zero-temperature pseudo-spin susceptibility 
profile Xzz (a;) = k--{x)/T and b) the Jacobian- weighted pro- 
files 'inx'^Xzz{x) in unit of 3A'^+/4eF+ is shown as a function 
of dimensionless position x for fixed interaction parameter 
l/{kF+as) = 3.0 at population imbalances P = 0, P = 0.2, 
and P = 0.6. 



and imbalanced Fermi-Fermi mixtures, we summarize 
our main results and state our conclusions. 



VII. SUMMARY AND CONCLUSIONS 

We derived general relations connecting the isother- 
mal compressibility and the spin susceptibility to fluc- 
tuations in the particle numbers for mixtures of equal 
mass fermions with population imbalance and for un- 
equal mass fermions. This derivation produced a gen- 
eralized fluctuation-dissipation theorem for Fermi-Fermi 
mixtures, that can also be applied to any other binary 
mixtures such as Bose-Fermi or Bose-Bose mixtures, and 
that is also valid beyond the local density approximation. 
The generalized fluctuation-dissipation theorem was used 
to connect theoretical calculations of the isothermal com- 
pressibility and spin susceptibility to experimentally ac- 
cessible density and density fluctuation profiles of bal- 
anced and imbalanced mixtures of equal mass fermions. 
We described the isothermal compressibility and spin sus- 
ceptibility for translationally invariant continuum sys- 
tems, as non-translationally invariant continuum systems 
in the presence of a trapping potential. Using the lo- 
cal density approximation, we obtained expressions relat- 
ing the local compressibility and local spin susceptibility 
to local fluctuations in particle numbers, and described 
their spatial profiles as a function of population imbal- 
ance. Lastly, we argued that discontinuities in these ther- 
modynamic quantities can be used to identify the phase 
boundaries between two qualitative different phases such 
as the superfluid and normal states. 
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Appendix A: INVERSE FLUCTUATION 
PROPAGATOR 

In this appendix, we present the elements of the inverse 
fluctuation propagator matrix F~^(q, where a;„ = 
2'KnT is the bosonic Matsubara frequency at temperature 
T . The diagonal matrix element reads 



Fn' = ^ + E |rkPC/q/2,k + T^, 



nintcr 
l/2,k 



and the off-diagonal matrix element 

k 

where C/q k and Vq k are the coherence factors 



1 - 



q-k,+ 



q,k 



-E-q+k / V £^q-k 
_ |Aq+k||Aq_k| 



4-Eq+k£'q-k 

and the intra-band and inter-band transfers 

T-iintra All , a22 

^q,k — ^q/2,k + ^q/2,k; 



rpintcr A 12 i a21 

^q,k — ^q/2,k + ^q/2,ki 



with 



"^q.k 



q— k.i^ j 



np{£, 



-'q— k 



1^ 



q+k,!/' 



(Al) 
(A2) 

(A3) 
(A4) 

(A5) 
(A6) 

(A7) 



Appendix B: Compressibility Matrix Elements 

In this appendix, we present the explicit expression of 
the pseudo-compressibility matrix element Kap , which is 
a partial derivative of Na with respect to fip with fixed 
T. Since |Ao| depends on implicitly, we need to add 
an implicit partial derivative of |Ao| with respect to fi p. 



Ka£ _ / dNa\ 



aiAoi^ 



d\Ao\' 



(Bl) 



But, since the uniform supcrfluid phase satisfy the or- 
der parameter equation (9r2/9|AoP) = and Np = 
— {dfl/d^p) simultaneously. Taking a partial derivative 
of order parameter equation with respect to ^p, we have 



d 



a|A^ 



9|Aop 



0. 



The first term can be rewritten as — (9iV^/9|Aop), lead- 
ing to 



g|Ao|- 



dN, 
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d\Ao\ 



Thus, the implicit partial derivative of |Ao| in Eq. (jBip 
can be expressed in terms of explicit partial derivatives of 
, Np , and the thermodynamic potential fi with respect 
to /ia, /i^, and |Ao| as 



T - [dtipj 



dNa 



9|Ao|V V^IAo 



_dH)_ 

9|A^ 



JB2) 

From the expression above, we notice that kap = Upa. 

Once we have solutions of the order parameter equa- 
tion together with the conservation of the particle num- 
ber of each species for given values of l/{kp+as) and P, 
we can evaluate the pseudo-compressibility matrix ele- 
ments kap as a function of l/{kp^as) and P. For in- 
stance, in the case of s-wave superfluid, the explicit ex- 
pressions for Eq. (|B2p are given by 



dNa\ ^ 



aiAoP 



E 



Aap{T) 



Ba{T) 

2E 



E 



k 



CjT) 
4E^ 



(B3) 



(B4) 



(B5) 



Here, we used the same notations X^ = 71^(61) ~np{£2) 
as in Eq. (|52p . And the second terms in the above equa- 
tions are given by 



v=l,2 



Ba{T) = 



C{T) = E 



\ d£. 



np{£u) 
d£. 



(B6) 



(B7) 



(B8) 



where = \ (7a (-1)" + i+/E) with 7t = 1 and 7^ = 
-1. 

At zero temperature, the Fermi distribution can be re- 
placed by the step function as np{£^) 9{—£,y), leading 
to 



X-\T=o^0{-£i)-e{-£2), 



(B9) 



and the corresponding the partial derivatives are by the 
Dirac delta function as {dnp{£,^) / d£„) — s- ~5{£y). For 
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convenience, let us introduce a dimensionless variable z = 
{k/kp,+Y, fl± = fJ,±/eF,+ , and Aq = Ao/eK+, leading to 
— Ei,/tF,+ ; where 

8,{z) = m{z - + (-l)-- ^{z- fl+y + |AoP. (BIO) 

Then, (5(15^) — (5(£^(z))/eF.+- As a next step in the cal- 
culate, we use the relation 



S{z-zi) 



(Bll) 



where £l{z) = m + {-ir{z - ^ {z ^ fl+f + \Ao\^ , 
and are the solutions of £^{z) ~ given by 



zi 



fl+ -mjl_ + (-1)'/D 



(B12) 



where D ~ [jl^m — — (1 — m^)|Aop, and m = 

(1 — mr)/(l + TO,.) with TOr = m-^/m^ the mass ratio of 
the spin-up to spin-down particles. 

Now, we can evaluate Aq,^(0), 5^(0), and C(0) by con- 
sidering the signs of D and z;. For instance, if £> < 0, 
A„/3(0) = Ba{Q) = C(0) = 0. But, if £> > and zi > 0, 
then we have 



3A^+ 
Aep- 



J2V^i\aMzi)WAzi)], (B13) 



^ B^jO) _ 3N+ ^ 



bgjzi) 

2E{z,) 



W.izi) , (B14) 



1. BCS limit 



In the BCS regime, it is known that Aq = Aq/ef 
(8/e2)e-('^/2)l^l and fi/ep ^ 1 with A = l/ikpa) - 
— oo [23, leading to (A//.i) <C 1, and 



9|Aop _ 2|Ao| 
9/i 



(C2) 



Now the susceptibility x — (dN/dfi) can be evaluated 
using Eq. (|C2|) . 



97V _ /5iV 

= E 



9|Ao| 



1 + a/A* 



El 



|Ao| 



(C3) 



Here, we introduce a dimensionless variable x = 
{ep / ^J^ik/kpy^ 1 leading in thermodynamic limit to the 
summation over momentum 



E = ^ 



(pk 

(27r)3 



C / dx vx, 
Jo 



(C4) 



where C = (3iV/4) (^/ep)^/^ ^ (3Ar/4). Then, we have 



dN 3iV , r ,2 /""^ , 

a;7 = 57l^»l i„ 



(i-l)2 + |A„P 



3/2 



(C5) 



By changing the variable a; to a; + 1 and separating the 
integration range in two parts 



dx 



dx 



x + 1 



3/2 



, (C6) 



E 



g(0) 



E 



AE^zi) 



W,{zi) . (B15) 



where a^pizi) = U'^{zi)U^{zi), ba{zi) = U^{zi), and 
W,{zi) = l/mzi)\. 



Appendix C: Analytical Result of {dN/d^)T for the 
Equal Mass Case at P = and T = 

In this appendix, we present analytical results of 
isothermal susceptibility {dN/dn)T for equal mass and 
population case at T = in both BCS and BEC regimes. 
The total number of particles = -I- A^^j, at T = is 



^ = E 



Ek 



(Cl) 



where = k'^/2m- ^ and E]^ = y/^l + |Akp. For 
s-wave superfluid, Ak — Aq. 



we can perform this integration. The integrand in 
the first region becomes l/{x^ + \Ao\'^)"^^^, in the 
|Ao| < 1 limit. Noticing that /^^ dx {x'^ + d2)-3/2 ^ 
{2/d'^)/VTTd?, the first integral is 

2 



|Ao| 



1 



AoP 



(C7) 



The integrand in the second region can be approximated 
by (x -I- l)3/2/2;3, since |Ao| ^ 1, leading to 



I2 = 



dx 



{X + 1)3/2 ^ 1 



7\/2 + 3sinh"i(l)] (C8) 



Substituting Eq. i^J} and dGSj) into Eq. (|C6l), we have 
(dN/d^i) in the BCS regime 



dN _ 37V 
9^7 ~ 27^ 



-7r/(fcF|a|) 



(C9) 



where 7 = [7\/2 - 4 + 3sinh"^(l)]8/e4 ~ 1.2519. 
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2. BEC limit 

Next, we consider the BEC regime, where kpUg — > 
0"'". It is known [l^ that the order parameter 
|Ao| — (16/37r)^/^eF /^/kpol and chemical potential /i = 
—E},/2 + (2/37r)ei?(fc_Fas), with E'f, = l/{ma1) a binding- 
energy, leading to 



I A, 



,2_ 4 4 



and 



2 37r 2m' 



(CIO) 



(Cll) 



For a given scattering length , we can perform the par- 
tial derivative of |Aop with respect to /i, by parametriz- 
ing the chemical via A:3 ^ N/V, leading to 



a|Ac 



9/i 



9|Ao| 



dkl ) \dk^ 



(C12) 



Taking A ~ \/{kpas) S> 1, we can rewrite |Ao| and the 
corresponding derivative as 



(C13) 



Introducing the dimcnsionless variable x = {k/kp)"^ /\^ , 
in the thermodynamic limit, the momentum sum reads 



(CM) 



Notice that the chemical potential in this large A limit is 
of the order of A^, since Eb/2 = epX^, and 

a ^ epX^{x + l), 



= epX^.y{x + l)^ + {m/3TT)/X^^ + O(A-i). 

Next, we consider the order of magnitude of the explicit 
and implicit derivatives. First, we notice that the explicit 
derivative of |Ao| with respect to /i is of the order of A~^ 



dN\ _^ |Ao|- 



0(A-2). 



(C15) 



Second, we establish that the implicit derivative is of the 
order of A 



0{X). (C16) 



Thus, in the BEC regime A ^ 1, dN/d^ can be approx- 
imated by the leading order term of the implicit deriva- 
tive: 

f .^±Ard.V^(. + l)- (C17) 
Noticing that dx ^/x (2- + l)^^=7r/2, wc have 



dN iN TT 



9/i 2€p kpa 



(C18) 



Since the total number of particles is = + iV^ = 
A:p/(37r^) and the Fermi energy is ep = kp/2m, we can 
express the derivative 



1 
V 



ON' 
9f^ / T 



(C19) 



in terms of the mass and the scattering length of the 
fermions. Since in the BEC limit all fermions are paired 
into tightly bound molecules a direct comparison to the 
results of weakly interaction Bose gas is possible. For a 
weakly interacting bosons, it is known |26j that 



1 fdNE 



Air \ as J 



(C20) 



in their superfluid state at T = 0. Given that the num- 
ber of bosons is Nb = N/2, and that {dfj.B)/{dfi) — 2, 
a substitution of the boson mass ms = 2m and its ef- 
fective scattering length as ~ 2as into Eq. (|C19|) . leads 
to {m/a) = {niB/aB)- Thus, the density susceptibil- 
ity for a Fermi superfluid in the BEC limit recovers the 
the expected results for weakly interacting Bose-Einstein 
condensates and lead to the relation: 



dNB 

d^B / T 



9fJ. /T 



(C21) 
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